{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Efficient handling of invalid simulation outputs\n",
    "\n",
    "For many simulators, the output of the simulator can be ill-defined or it can have non-sensical values. For example, in neuroscience models, if a specific parameter set does not produce a spike, features such as the spike shape can not be computed. When using `sbi`, such simulations that have `NaN` or `inf` in their output are discarded during neural network training. This can lead to inefficetive use of simulation budget: we carry out many simulations, but a potentially large fraction of them is discarded. \n",
    "\n",
    "In this tutorial, we show how we can use `sbi` to learn regions in parameter space that produce `valid` simulation outputs, and thereby improve the sampling efficiency. The key idea of the method is to use a classifier to distinguish parameters that lead to `valid` simulations from regions that lead to `invalid` simulations. After we have obtained the region in parameter space that produes `valid` simulation outputs, we train the deep neural density estimator used in `SNPE`. The method was originally proposed in [Lueckmann, Goncalves et al. 2017](https://arxiv.org/abs/1711.01861) and later used in [Deistler et al. 2021](https://www.biorxiv.org/content/10.1101/2021.07.30.454484v3.abstract)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Main syntax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sbi.inference import SNPE\n",
    "from sbi.utils import RestrictionEstimator\n",
    "\n",
    "restriction_estimator = RestrictionEstimator(prior=prior)\n",
    "proposals = [prior]\n",
    "\n",
    "for r in range(num_rounds):\n",
    "    theta, x = simulate_for_sbi(simulator, proposals[-1], 1000)\n",
    "    restriction_estimator.append_simulations(theta, x)\n",
    "    if (\n",
    "        r < num_rounds - 1\n",
    "    ):  # training not needed in last round because classifier will not be used anymore.\n",
    "        classifier = restriction_estimator.train()\n",
    "    proposals.append(restriction_estimator.restrict_prior())\n",
    "\n",
    "all_theta, all_x, _ = restriction_estimator.get_simulations()\n",
    "\n",
    "inference = SNPE(prior=prior)\n",
    "density_estimator = inference.append_simulations(all_theta, all_x).train()\n",
    "posterior = inference.build_posterior()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further explanation in a toy example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sbi.inference import SNPE, simulate_for_sbi\n",
    "from sbi.utils import RestrictionEstimator, BoxUniform\n",
    "from sbi.analysis import pairplot\n",
    "import torch\n",
    "\n",
    "_ = torch.manual_seed(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will define a simulator with two parameters and two simulation outputs. The simulator produces `NaN` whenever the first parameter is below `0.0`. If it is above `0.0` the simulator simply perturbs the parameter set with Gaussian noise:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulator(theta):\n",
    "    perturbed_theta = theta + 0.5 * torch.randn(2)\n",
    "    perturbed_theta[theta[:, 0] < 0.0] = torch.as_tensor([float(\"nan\"), float(\"nan\")])\n",
    "    return perturbed_theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The prior is a uniform distribution in [-2, 2]:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior = BoxUniform(-2 * torch.ones(2), 2 * torch.ones(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then begin by drawing samples from the prior and simulating them. Looking at the simulation outputs, half of them contain `NaN`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6ee9400b22f4437ca720a451d7af7bd0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Running 1000 simulations.:   0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulation outputs:  tensor([[ 0.0411, -0.5656],\n",
      "        [ 0.0096, -1.0841],\n",
      "        [ 1.2937,  0.9448],\n",
      "        ...,\n",
      "        [    nan,     nan],\n",
      "        [    nan,     nan],\n",
      "        [ 2.7940,  0.6461]])\n"
     ]
    }
   ],
   "source": [
    "theta, x = simulate_for_sbi(simulator, prior, 1000)\n",
    "print(\"Simulation outputs: \", x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulations that contain `NaN` are wasted, and we want to learn to \"restrict\" the prior such that it produces only `valid` simulation outputs. To do so, we set up the `RestrictionEstimator`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "restriction_estimator = RestrictionEstimator(prior=prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `RestrictionEstimator` trains a classifier to distinguish parameters that lead to `valid` simulation outputs from parameters that lead to `invalid` simulation outputs "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training neural network. Epochs trained:  35\r"
     ]
    }
   ],
   "source": [
    "restriction_estimator.append_simulations(theta, x)\n",
    "classifier = restriction_estimator.train()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can inspect the `restricted_prior`, i.e. the parameters that the classifier believes will lead to `valid` simulation outputs, with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The classifier rejected 51.6% of all samples. You will get a speed-up of 106.5%.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAJgCAYAAABr3Xx6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtvElEQVR4nO3df4yd51kn/OtObSfuTJKxcf0jdkPabNKo0MZrNu3upgkEBdTwM2yJlFcFVEQF7Asv21VhNwvalkViyb5Sl92l3S1QaLba0laFErVq3ZKqWSVtoHLxmxAIxYmbKPXvdW3HHtvN2Ph+/5gxmvq6Jz62Z+bMzP35SFbsO+c59/0858zMd55zPddTaq0BANCzy4a9AACAYROIAIDuCUQAQPcEIgCgewIRANA9gQgA6N6yi9jGdfrQnzLsBcyWb//A/5u+h43sXJ4et2rH36exwze+LI1t/P7n09iRB16ZxiaubB/C8TedSGPrP35587Fp2w15PS1HbziTxpaND/aSvnxvftypKwfaNI5ff6o5vmJ//tHzjz64N43tefM1A81zZHOeZ9loHlu9dWXe9sb8fGu35+PVOtbrvvxCGnv27qvS2PIZjvXInvzj9Pg1+bHLjzU3T1qvy93/16Np7MGP3Ja3Hc1raR2HPbcP9r655pF2VDiwJZ+HmVh3Oo2NPZ6/Jld/9cU0duim/LXSetzDn7/vvAt3hggA6J5ABAB0TyACALpXLuLWHWqIoD9Lpobo9//utvQ97L2/8y/S405syN/qWrUOlx3PtSVrtufDteqpo831nLr6ijR2xXPfSGMHb92QxlaMD1bjcdm6b6ax0S++PI0NWs80MZp/l27Vvrx4c36+mZ6zVVsytiNv2zqOrbqdTQ9PpLFdd6xIY6cbtTNXX384jR3el+do1Z61zFRL1aqTab0GrdqnVk3aimN5X1o1Ui2tWrHW10CrHqpVC3UhWu/jVs1Wa/9a77t129QQAQBcFIEIAOieQAQAdE8gAgC6p6gaGMSSKaq+48770/ew59+emzBe/kQuOm41yBu00LPV6DGiXUw8aOFqq2C2VUzcal7XatbYKiY+/sTqNNYqwG01R5yp6HjQJoCtRpitAuW1j+ZGj63XoHUMWwW93/3ux9LYxz/1pjTWep1aheStAvaIdiPF1nts0GL31nvpzPUnB1pP6308aNF+y0yNGVsF+c3HNYrGW8eh9XXa8tRv/WtF1QAA5yMQAQDdu5h7mcElufX+L8TuIydj49jK+NJ93zvs5QCAM0TMv91HTsZz9/9g7D6SP9u+ELfe/4W47r5Px633f2GWVgZAr5whYtE6G6yuu+/Tw14Ki8j+W3Ix6mU7cwFoq/D3xZtziH/1uoNp7Mhj+W73K17IXa4jIk6P5sLjVofmVgF1q7i51Y15+XjetnW3+1X/ZTRvuzFve+C2XNj8mt/NBa873zlYJ+eIiMPjuTi29Rq0OoO371ifi9jb3bDz2Kd/L98NPhpdm1uF8iu25fVNXD1DQf2G1jmJ/LqcHs/HsVWc3LxrfKOAulVQf2Ykr7FVrL5y32BF1RON1y6iXSzd0lrjNc1C8tm7zssZIgCge84QMTQbx1b+w9kd9UQADJNAxNBMD0A+9gJgmHxkBgB0zxkiFgWX6jNb1m17MY21CpFbnZPH9+YC1a+9aU0au3bAbtEREWdGcoHywS25MHflvvz767ovv5DGDr/2qjR25MbBCk9bx6FVuL1sNK/5wLtzQe/pfe0uwu3O0oOssN0Bec/t+XHjG/IxHPRGC6N7c4Hx0RvycWgVaa/dnvetvZaITQ/n92LrfdIqWH/6J3Kb6xuv35vGdsSGvMbG8V8xnvdlfMNgnbjP7L8ijY3saRdPt7pNt4rGW2v87nd/KY195OFb01jra2UQzhCx4LQup5+tS/UBoMUZIhacYV9Of/ZsVIRib4BeCERwjrOBLEKxN0AvfGQGAHTPGSKWlKVYfL0U92mYDr8jt21e/qVvS2OtQthTuY41Vm9dmZ/vhaNpbMWxXLAcETGyMxeUnlyfC7qXH8vbtgprr3o6F7NOrMtF0GOP53lb3YFXNOZtFcEeboy15oiImGgcx9b+nVyfi6BbHZoHLZZuOXRXrktsFQm3bBqweP6az+5pbn/w1lzw3HLq6ryeVmfp8m+uTmMjd+bXYNVTuRh/151529Zrcu0H8tfFoZvyazJTkfz6RrfpPbfn9+ypew6lsQc/kjuIr9mTX/s1X2oc799qr2c6Z4hYUuaj+Pps0fd83UdNQTnA3HOGCC6QGiOApUcgglm20D/ichUdQOYjM5hlC/0jrrPrW8hrBJhvzhDBkCz0M0lL1StGjqex50ZXp7Hlx3KhZ6vYed22XLDcKoI9uKVd+Nsqgl61Iz/2QGP7Ndvztj/7q59IY+/9nX+R19gobL76+sNp7Pi+XHC+aWv+XfrwjbnY9sSG9j63ul+3HjvoelrHsFUg3jpeqx7Mr19EHjt+7Wga239LLhBetSMXOz/z0+3i6dOjrX3OxcTPvykX7o99Mc/dKoxuvWd3vjMXWp/ZP9h7e3xjvjig9dq19i0iYmRPfu+s2Z4fdzBWpbF8FNpfVwe3rG/OfT7OEMGQnFuLNF9F2gBkzhDBkE0/O6RIG2A4BCIWhI1jK/8hDGwca50YXZjOrtvHXgCLm0DEnBq0Tqb1/2YrJM1lrc7Z51uoZ3ZcUQYwGIGIOXUpN2od9If3+YLTsG8WO925AWWuna9nkjNcAJMEIha92fpBPv1M0lnnhq1LnWt6QFkIFvoZrrlw5IFXprGxxuPG3vZ8GjvxWN62daVR68qemaz+6otp7Pm35yuVRr/48jR26p5vpLFBryhbty3Pu6txtd2Zm0+ksfFjeS23/NiTaWzbn74uTxztq5LGduTHHYm8nrWNq7haTl2Zr3pr3eKjdWXWi419jsjH6/In8nG4ENc8ktdzuHEV3WjjFhqje/NxGI+8zycbF1y1bk2ycl++xqp1G5K2vB+t54to39Jj4/fnr7WDf7VpoJlbVw62bvsxCIEIprTCioJngD4IRMyal6rVma16oMVafD2oi/kISz8jgEsnEDFrXqpWZ7Z+UC/1H/gX8xHW2eN+9qazEQqoAS6UQASLxPnqmXy8B3DxBCJm3VL/WGtYBJ7ZcfyaRhFm4zYDh/avSWOtMtFNn38hje1/Yy7UXfuWXDgaERGvz0PX/ua6NPa1e/NtFJaN56+vfIOJiOPX5233Ry6YXT6et1318fy4icbxeuY/vTZvG+0C6ANb8pFsvS5rt+ftW9u2CrJbDt2V793XKlbfuO5gnvdPrk1jrduDbPz+3Wns9AwFwq3bnTSL3e/It8uIRgF1q1j9uk/l1/65H8637lhxay7QX/7xXNS+Yjzv89q35H3e/Wf5eEW0v14O78kXK+SvvogjN7ZvB3Kum9fuGehx5xKImHU+qgFgsRGI6Np89wUCYGESiOjafPUFavU4AmDhEIhgHiy0howAfCuBCBaghVCYPv3jxKUU5lYN2Ol417pceBrrTqehne/Mj7tsZy7+PPHejc15JkZzkfCRO/Lj1j7a2PbKPHersHbF/vytvlVIPqgjN+ax06O5KPqqp9vdijc9nAuHW12RxzfkwuFW4XerU3Jr/1qdpVuF0Ssa3cyj0e37hv+VW0g/e3dj2xmO9fJGB+qJq/Nrdbqx/akr8/E+c30uGt8V+fvHyn1524lGh+yJRlfp441i7gONDu6nZ+jW3rrgoKX1urQ6e7feI49/qNEh/Q3nn1MgggXoUgvTpzd4vFDTP95bSkEI4KUIRDALWmd0hnmW51IClY/3gB4JRDALWgFE+wGAxaP9AS8AQEecIQK6suuuXKx52fFcmLn20UG/PebHtYqOV7yQC7IjIvbcngujl43notdWJ+eTjcLVVuHpntvztq05rvvnX09jB/bmjsNjO/Icx6/Jv1+3CrwjIo7ekPd55b7mQ5NWR+WxkeNpbLTR7bvVoXnlvlbX7Dxvq4P03/1cLtJ+9Ucn0li703TE+JtOpLEVWwcrgn7x5rxtq+v2qUYxeEurCL1VzH1mJF+UsGlrPoatLtwRESuOtd47gxbk57lbx3B14xgOQiBiSRqkpsdHWgCcJRCxJJ2vpqfHe4FdypVnAEudQAQDWuyBwhkxgJkJRDAggQJg6RKIuCRujspi0+ravHZ7Lk6eaHRePnBbLowe2dnoaN3onNzqxDwpz90qZp1odMluzX1gS962VZTbKhI+siN3HD5x42BdkltFvqfHW8em/Rq0uja3HH9idRqre3OX5RN3DHYcBi0mPvyO3CJ70wfH0tihm/I+v3xvGppcY6OLdOs91jpeLeXNueB803/I2z5791UDPV+rQH/k+fw6v+xwLmrfddcrms95cn0utj4zkvd57PF8HFsXB1zVKCQ/dFde4yAEIi6JJn4ALAX6ENGlxV4PBMDscoaILqkHAmA6Z4gAgO45QwR0pVkw2+iq2yq2bRUxb/z+5/McjS7JMxndO9i34T235zW2CpGXH8uFp61i20PHciHy0RtygXerk/N7f/b9aexfv+fn82Jm0OoufGRdPratDsinrszHYd2XX0hjT/9EbtHcKp6fqaPyuepn8/Ea35AfN7o3d1OeqUv5imO5g/WK7Xmf991zMs/TKCY+viGXADz9E/l9fNXTjeLyRkfr4z99JI2NNIq0D96aD8TY43mOmRzZnI/Zkc2n0tirP5pfv/235IsVWscm7j3/OpwhAgC6JxABAN0TiACA7qkhApo3wwXoiUAEi9BsB5ie2hBc96lcrLnrjlzc2uoM3fqW+cxfbUpjK2/JJ99bnaEjIvbd0x5P9l+RhlrF0i/sXJXGrmoUBF/z2T1p7OgN69NYq7j85/7o59LYPT/7aBp76D1vSmMREZc/kYteX90oTn/m9nxsz4zk1+/4jtE0tmw8F/W+7T8+mMZ+7z/+WHON52p1P371uoNp7MgDudv3gS35/RUR8fK9eY3Hr8lF3td+II89//a8nlYxcaur+Ipj+TWNyI8baXTi3nXnYEXoJ9fnAuiIdvfrdV/O+3L4tbmb9rJf+3oaG2kc7xXj7bnPRyCCRainAAMwH9QQAQDdE4gAgO4JRABA90qtreKql3TBG7B0XXffp93tvg+Dt51d4F71X9+Tvoet2d4qbs1jrULRQbdtFSdHtAtrT2zIjx3bkbdtzfPizblA9bKdufC+9XwtE42i3FZX49axaRXQRkQc2JJ/F291EL/q6fy4Vjft1jyt7tATV+ey2UG7lLe6XLf2o1WMv/bRdrnuqqeOprFn787FxK3nHHs8d/ZuvW+WN4rLW/s36PurdXFAq1v0qh25+/RMBn0NWl8rrQLxQ3flzt477/21834Pc4YIAOieQAQAdE8gAgC6JxABAN3TmBHoSqtQ98Btufvxiv3522Ory/X4xtyFePmxPO+p3Ew5ItqFsJdi/cdzgev4hlxPenDLYAXirQLqltZxHd/QfuzIzbnD9rHxRsf1p3Pn5ZX7Gq/flnwMJ9a1ft9vdTDOY6/+aB5rdTN/1YO5KHrXnVensdHd7W7krQLqTQ9PpLFDN+XXdPVX83Oe2NDuiD2Ig41j2LqGqvV+bxVQz9SBffXW/Drf8mNPprHHP/S6NNYqoD51z6E0NtrozB73NpfzLZwhAgC6JxABAN0TiACA7glEAED3FFUD3Wt1Em51aG515G11hm4Z/WIuEI6IuOaze9LYU/e9Io0dH8+diVtaxazXvycXgx+9IVdLD9qdu1Vc3io6PnN97hgcEXHZE6vT2OWNjsqtgvORm3MR7cmdq9LYqvW54Lk2im1bReOHbspjrQLj/W/MBdStTs6touiIdnH5c6N5X5aN57nLm8fT2PIv5f07fn1+rUZ25vdSqwt069gcvyaPrfn4/jS26678Ho6IGN2di8af/9Ub09iRewe70GHtB8fS2EzF/OfjDBEA0D2BCADonkAEAHRPIAIAuqeoGujKkc25WLNVVL32Lc+nsd1/dm0au/yJXCzdKkQ+ekOrS3LEkc3t4tNztboBH9iSf6e99gMvS2O77sxrPDOSj8PyY/n5Vu1odIG+utHF+59/PY3t2Nmubh1pFFC3jlmrc/bBRtFxq0v2xL5cYDzS6nTcKJ4f9D1y5Mb8fEdvyAXLraLoiIiXtzoqN94nrYLnurdRIN4oQh97PK+nVax+/PrTA23b6pD91f9nfRpbsT+vOSJi/y2N4u3RvJ4VuU676fCN+f3eKiQfhDNEAED3BCIAoHsCEQDQPYEIAOieomqgK63i2Nf+wl+nsScO5Ja8I3ty8eehu3I35jVbV6axg1vahbWt7rvLG0XH++7JHbFXN+Z5/u25+Hr9x/PYyfW5GLVVTBwxWFHu2AOvTGMjjc7XEe1uzuMbc6frVrfwZblBc1OrsHb5sbwvrU7jl+2/YqC1LG+spfXaLT/WXuOKRpH3Df8rP3jnOxvFzY3XPiLPvfmnnkxjjz30nXktjfdh6/3QKhpvFb+fuid3FI+IODae132mcbxX7svna1rF162vydbrPAhniACA7glEAED3BCIAoHsCEQDQPYEIAOieq8yArhxvXPm07U9fl8baVwblK1pGv5hvi3FwS/s2HS2nG1fOrN2etz8ceZ4V4/nqsStH81Vve27PV/Esa1wNFeOD/Uho3Upi4sr8uNbtOCIiDt10eZ76Tflqr9a+tG530brqbdlovkJq9Vfzek5syFc9Xd44Nq1br5wZycd/ZOdgt7uIiHjuh/Njj19zdRq7cvQbaezIja2rzPJxeOp9+YqylzeumCtvznOs/eBYGpsYzdu2vqZOPbG6sb72sV1xa5778MhVaWzT1nwOZ8/t+fmuerp9deP5OEMEAHRPIAIAuicQAQDdE4gAgO4pqga6smpHLoTdd08uej21c7BbI7SKUa9qFP6ufcvzzfWUf5OLaD/7qQ+nse/8r/93Gjt8Y779xi1r96Sxvf8hf6s/8O7T+fn25ULWtdsH+71511256Lh1m5SIdgH2ZY3jvfbBXBi968687aaHJ/LjIj/f/ltaxba5ELlV5Hu8cWxWrT+axka2jjXmzUXkERGXrcuF5KseyY89MJoLlFuF7acax/XUPXlfTn0pvz/zjVPa76/W10/sye+RIzc2njDaFyu8YuR4Gnvh+Ko01irobhVQz1TEfj7OEAEA3ROIAIDuCUQAQPcEIgCge4qqga60OtvG/tzJeX2jW3SrcDiXfraLW3f/2bXN9bz4zlxYe/3Hfj6NXd3q5tso9H3sodyZeOTdh9LYsfFcdNzqsrzvnry+VgH0ZcfzcZ2psHZsR6OQuVFs++zdjW7FD+eC2VbH56uvz/s80SgmbnUF33d93r+xx/McNfLzReSi43bX83bh/p7b87FZua/R3flYo2v63laH87zGVmH0ruvzsX71tnysv3ZvPo/ymt8dT2OtjtsREUc250L5tY3HtTqpH9zS6Eiep47VX21OfV7OEAEA3ROIAIDuCUQAQPcEIgCge4qqga6s3Jd/DxzZk4s199yet20VHS9/NHcRPnxbLhxtbRsRcfkTL09jp0bzel7Ymcu3W7/RtoqEx/fmwtrLNuQ5Tq7P2171xby+Vnfu40/k47C8URgbEXGgcXxW7B/sx9H4xtxTuVV0PLEv7/Omz7+Qxva/MRf/rt7a6lKej9fElXneVnfn1nGNiFizPW8/ujsfm1bR+Mn1+flW7stzH78+P9+p0XysW4XRrWMz9niet1X8PrEuzxsRcdnxvMbWBQcTjXW3rH0kfxUcuqndGfx8nCECALonEAEA3ROIAIDuCUQAQPcUVQNdaXUN/k+//rtp7PPHviONPfSeNw00x6at+XfNXXe1i0TXPpq/DY+97etp7MCf5MLTVkfskedz4XDEaBo5ekMu6D0zkjsYn7oy78vJRoH3ykYB9bpGp+OIiNiWhw7dlLc/0Sj8bnUrvuaRXLTc6kjeKhJuzREx2FpahdYzFZK3jO6eSGOtovEzI6fT2LLRxvtpXy6Ab72/Vozn43Xq6tytvfX+ahVpt97vI8/nDucREU//RH7S1kUNI3vyulc9dbT5nOeaGM1F3oNwhggA6J5ABAB0TyACALonEAEA3Su1tgrFXtIFb8DSdd19n47n7v/BYS+DuTd4pegCd9c/+pWBvof9/aqRNHb82lyc3DIxmn/XbBWyzvTYllYB7q47cgHuqx7Mhaetgtn9t+Ruvq3i1pZWcWurYPnoDYN3aG4VS7eKvFudjq96Oh/DVhF0q+C51b26dbxaLv//vpbG/v7w4YG2jYh42apcnH4h2zO4h858/Lzfw5whAgC6JxABAN0TiACA7mnMyEW59f4vxO4jJ2PjWOuu0ACwuCiq5qIopu7Okimq/r7L7vE9DDqjqBoAYAACEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQvWXDXgCLx633fyF2HzkZEREbx1YOeTUAMHsEIga2+8jJeO7+Hxz2MgBg1vnIDADonkAEAHRPIAIAuicQAQDdE4gAgO4JRABA9wQiAKB7AhEA0D2BCADonkAEAHRPIAIAuicQAQDdE4gAgO4JRABA9wQiAKB7AhEA0D2BCADonkAEAHRPIAIAuicQAQDdE4gAgO4JRABA9wQiAKB7AhEA0D2BCADonkAEAHRPIAIAuicQAQDdE4gAgO4JRABA9wQiAKB7AhEA0D2BCADonkAEAHRPIAIAuicQAQDdE4gAgO4JRABA9wQiAKB7AhEA0D2BCADonkAEAHRPIAIAuicQAQDdE4gAgO4JRABA9wQiAKB7AhEA0D2BCADonkAEAHRPIAIAuicQAQDdE4gAgO4JRABA9wQiAKB7AhEA0D2BCADonkAEAHRPIAIAuicQAQDdE4gAgO4JRABA9wQiAKB7AhEA0D2BCADonkAEAHSv1FqHvQYAgKFyhggA6J5ABAB0TyACALonEAEA3Vt2oRuUUv46Ir45B2sZ1JqIONjp/D3v+7Dn73nfIyKuqLV+5xDnB5hTFxyIIuKbtdZ/MusrGVAp5Su9zt/zvg97/p73/ez8w5obYD74yAwA6J5ABAB072IC0e/N+irMvxjm7n3+nvd9IcwPMKd0qgYAuucjMwCgewIRANC9CwpEpZS3llL+qpTyZCnlsVLKzXO1sFLKTaWUPy+lvFhK+eWXeNwDpZRnSymPT/3ZPMz1zNJcpZTy30opz0wd7y0zPO5/l1L+btq+r53DNb15aq5nSin3DWueUsrbSin/Z9o+v32u1jI13x+WUg5M9d+aU+ebq5TyPaWUF6bt+7vmeD2vLKU8XEp5qpTyN6WUfzWX8wEM04X2IXo2Ir671nq4lHJXTBZavnH2lxUREYci4pci4u4BHvsrtdY/nqN1nHUh67lUd0XEDVN/3hgR/yNmPs5vrbXOaY+YUsrLIuJ9EfF9EbErIraVUj5Za31qSPN8rNb6i7M590t4ICLeGxEfWiBzPVpr/aF5WEtExOmIeGetdXsp5cqI+MtSykOz/boDLAQXdIao1vpYrfXw1D//IiI2zf6S/mGuA7XWbRFxaq7muBDzvJ4fjYgP1Ul/ERFjpZQN8zDvTN4QEc/UWr9Wa52IiI9OrXGxzjOwWusjMRmGl9Rcg6i17q21bp/6+7GI+NuI2DjcVQHMjUupIfqZiNg6Wwu5RL859dHSb5dSLh/2YmbBxoj4+rR/74qZfxB9cOrjk39fSikLYD3zMc9bpl7vPy6lvHIO1rGQ/bNSyhOllK2llO+Yr0lLKddFxD+OiC/P15wA8+miAlEp5Y6YDET/dnaXc1H+XUTcFBG3RMTqWBhrmi9vrbW+LiJum/rzk0Nez3z4VERcV2t9fUQ8FBH/c8jrmU/bI+Lba603R8TvRMSD8zFpKWU0Iv4kIt5Raz06H3MCzLfzBqJSyi9MK+K8ppTy+oj4QET8aK31G7O5mHPnGmSbqdP6tdb6YkR8MCY/dhnaemZjrojYGxHTz3xsiojd525Ta9099d9jEfFHMYv7fo7dg6xnPuaptX5j6rWOmHwfftccrGNBqrUerbWOT/39MxGxvJSyZi7nLKUsj8kw9OFa6yfmci6AYTpvIKq1vq/WurnWujkmi7A/ERE/WWvdMduLmT5XrXXPINucra2Z+rjo7oiYtauBLmY9szFXTP7m/1NTV5v904h4oda6d/rjSynLzv4wnPqh9UMxi/t+jm0RcUMp5VWllBURcW9EfHIY85xTS/UjMVnX0oVSyvqzH4uWUt4Qk1+/s/pLyTnzlYj4g4j421rrf56reQAWggu9yuxdEfFtEfHfp74vn56rO3CXUtZHxFci4qqIOFNKeUdEvLbWerSU8pmIePtUSPlwKeUVEVEi4vGI+Pn5Xs8cTPeZiPiBiHgmIk5ExE9PW8fjU6Hp8oj43FQYellEfD4ifn8O1hK11tOllF+MiM9NzfWHtda/ma95Sim/ERFfqbV+MiJ+qZTyIzF5BdShiHjbbK9julLKRyLieyJiTSllV0S8u9b6B/M1V0Qsj4iotb4/In48Iv5lKeV0RJyMiHvr3LaavzUmP4Z9curMZUTEr06dnQJYUty6AwDonk7VAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gIkopv15K+eWpv/9GKeXOS3iuebs7PADMFoGIb1FrfVet9fOX8BQPRMSbZ2k5ADAvBKJOlVJ+rZSyo5TyxYh4zbTxB0opPz719+dKKb81dUuRr5RStpRSPldK2VlKaTbAXGh3bAeAQVxop2qWgFLKd8XkbTE2x+R7YHtE/OUMD3++1rq5lPLbMXn259aIuCImbxPy/jlfLADMA4GoT7dFxJ/WWk9ERJRSXuq+ZGf/35MRMTp1I9ljpZQXSyljtdYjc7tUAJh7PjLjfM7eWf7MtL+f/bdADcCSIBD16ZGIuLuUsrKUcmVE/PCwFwQAwyQQdajWuj0iPhYRT0TE1ojYNlvPPXXH9j+PiNeUUnaVUn5mtp4bAOaKu90DAN1zhggA6J5ABAB0TyACALonEAEA3ROIAIDuCUQAQPcEIgCgewIRANA9gQgA6J5ABAB0TyACALonEAEA3ROIAIDuCUQAQPcEIgCge8uGvQCAeVaHvQBgXpVBHuQMEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdM+9zAC4ZLfe/4XYfeTkt4xtHFsZX7rve4e0IrgwAhEAl2z3kZPx3P0/+C1j19336SGtBi6cj8wAgO4JRABA9wQiAKB7AhEA0D1F1QBEhCvF6JtABEBEuFKMvvnIDADonkAEAHRPIAIAuicQAQDdE4gAgO4JRABA9wQiAKB7AhEA0D2BCADonk7VAA2D3sbC7S5gaRCIABoGvY3FbN/uQsCC4RCIABYQ9xO7dELlJMfhwghEACwpQuWkxXAcFlJoE4gAYJqF9EN6qVtIoU0gAoBpFtIPaeaPy+4BgO4JRABA9wQiAKB7aogAmBMbx1am2hvFyRdGgff8EYiArrR+wLRsHFs5D6sZnpl+0M6m1g9txckXRoH3/BGIgK60fsD0yHGAbyUQAQzJfJylWQzcN47phvVRq0AEMCTO0kwa1n3jWJiG9VGrQASwCF1ILZQzKLwUZ94mCUQAA5rpVP4wDHp2yRmU2bGUr5hz5m2SQAQwoMX4w6/1g/zsOINzxdzSJxABLGFzEeLmoxh8IZ2NWwwU6F86gQiACzIfxeCL8WzcMCnQv3Ru3QEAdM8ZIgDmjY/CWKgEIoAFbimFCB+FsVAJRAALnBDBYrMYexsJRACwiCyGM4aLsbeRQAQAi8hCPsuymAlEACx5MzWonOmx9EcgAmDJc1aF89GHCADonkAEAHRPIAIAuqeGCAA6sRgu2R8WgQhgHrgbOQuB4vKZCUQAs2ym38LdjXz2OePBbBGIAGaZ38LnzzCP9UxhzOs/++bjWAtEAHARWj+MF/rtKRar+TjWAhEAM/KRFL0QiACYkY9/5sZivBv8UicQAcA8W4x3g19oBr0/3aAXMwhEAMCiM9tn03SqBgC6JxABAN0TiACA7qkhAoA55LYtkxZ6CweBCADmUOuKspaFHhgu1UJvKSAQAcACsNADw1KnhggA6J4zRAAwS5b6x14tS6VGSiACgFmy1D/2min8DNoNeiETiACAbzHTbTGWSvhpEYgAgG+x1M90tSiqBgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOieQAQAdE8gAgC6JxABAN0TiACA7glEAED3BCIAoHsCEQDQPYEIAOheqbUOew0AAEPlDBEA0D2BCADonkAEAHRPIAIAurds2AsAmE+llL+OiG8Oafo1EXFwSHP3Pn/P+z7s+Ye971fUWr/zfA8SiIDefLPW+k+GMXEp5SvDmrv3+Xve92HPvxD2fZDH+cgMAOieQAQAdE8gAnrze53O3fv8Pe/7sOdfFPuuUzUA0D1niACA7glEQBdKKW8tpfxVKeXJUspjpZSb53Cum0opf15KebGU8ssv8bgHSinPllIen/qzeZjrmaW5Sinlv5VSnpk63ltmeNz/LqX83bR9XzuHa3rz1FzPlFLum6t5BpmrlPK2Usr/mbbfb5/j9fxhKeXAVLuJOXW+uUop31NKeWHavr9rjtfzylLKw6WUp0opf1NK+Vcv9XiX3QO9eDYivrvWeriUcldM1hW8cY7mOhQRvxQRdw/w2F+ptf7xHK3jrAtZz6W6KyJumPrzxoj4HzHzcX5rrXWgS6IvVinlZRHxvoj4vojYFRHbSimfrLU+NcS5PlZr/cXZnn8GD0TEeyPiQwtkrkdrrT80D2uJiDgdEe+stW4vpVwZEX9ZSnloptfeGSKgC7XWx2qth6f++RcRsWkO5zpQa90WEafmao4LMc/r+dGI+FCd9BcRMVZK2TAP887kDRHxTK31a7XWiYj46NQaF/tcA6m1PhKTgXhJzTWIWuveWuv2qb8fi4i/jYiNMz1eIAJ69DMRsXXYi5jym1MfLf12KeXyYS9mFmyMiK9P+/eumPmH0AenPjr596WUsgDWM19zvWXqNf/jUsor52gtC9U/K6U8UUrZWkr5jvmatJRyXUT844j48kyPEYiArpRS7ojJQPRvh72WiPh3EXFTRNwSEatjYaxpvry11vq6iLht6s9PDnk98+VTEXFdrfX1EfFQRPzPIa9nPm2PiG+vtd4cEb8TEQ/Ox6SllNGI+JOIeEet9ehMjxOIgCWrlPIL0wo4rymlvD4iPhARP1pr/cZczjXINlOn9Gut9cWI+GBMfuQytPXMxlwRsTcipp/12BQRu8/dpta6e+q/xyLij2IW9/0cuwdZz3zNVWv9xtTrHTH5XvyuOVrLglNrPVprHZ/6+2ciYnkpZc1czllKWR6TYejDtdZPvNRjBSJgyaq1vq/WurnWujkmLyL5RET8ZK11x1zOVWvdM8g2Z2trpj4uujsiZu1KoItZz2zMFZO/9f/U1NVm/zQiXqi17p3++FLKsrM/CKd+YP1QzOK+n2NbRNxQSnlVKWVFRNwbEZ8c1lzn1FP9SEzWtXShlLL+7EejpZQ3xGQGmdVfTM6Zr0TEH0TE39Za//P5Hu8qM6AX74qIb4uI/z71Pfn0XN1wspSyPiK+EhFXRcSZUso7IuK1tdajpZTPRMTbp0LKh0spr4iIEhGPR8TPz/d65mC6z0TED0TEMxFxIiJ+eto6Hp8KTZdHxOemwtDLIuLzEfH7c7CWqLWeLqX8YkR8bmquP6y1/s18zlVK+Y2I+Eqt9ZMR8UullB+JySugDkXE2+ZiLWeVUj4SEd8TEWtKKbsi4t211j+Yr7kiYnlERK31/RHx4xHxL0sppyPiZETcW+e2O/StMflR7JNTZy8jIn516uxUXr9O1QBA73xkBgB0TyACALonEAEA3ROIAIDuCUQAQPcEIgAWrVLKr5dSfnnq779RSrnzIp/ngu6MztKjDxEAS0Kt9V2XsPkF3RmdpccZIgAWlVLKr5VSdpRSvhgRr5k2/kAp5cen/v5cKeW3pm4p8pVSypZSyudKKTtLKakB5oXeGZ2lxxkiABaNUsp3xeQtMTbH5M+w7RHxlzM8/Pla6+ZSym9HxAMx2bn4ipi8Tcj7X2KO6+I8d0Zn6RGIAFhMbouIP621noiIKKW81H3Jzv6/JyNidOrMz7FSyoullLFa65FzNxj0zugsPT4yA2CpOntX+TPT/n723+mEwIXcGZ2lRyACYDF5JCLuLqWsnCp+/uHZeNILvTM6S49ABMCiMVX4/LGIeCIitkbEtll66rN3Rv/eqULsx0spPzBLz80i4G73AED3nCECALonEAEA3ROIAIDuCUQAQPcEIgCgewIRANA9gQgA6J5ABAB07/8Hi+TTkSqTbx8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "restricted_prior = restriction_estimator.restrict_prior()\n",
    "samples = restricted_prior.sample((10_000,))\n",
    "_ = pairplot(samples, limits=[[-2, 2], [-2, 2]], fig_size=(4, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, parameter sets sampled from the `restricted_prior` always have a first parameter larger than `0.0`. These are the ones that produce `valid` simulation outputs (see our definition of the simulator above). We can then use the `restricted_prior` to generate more simulations. Almost all of them will have `valid` simulation outputs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The classifier rejected 50.9% of all samples. You will get a speed-up of 103.6%.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1f2d11a367bd4b89b301213ebcbf478f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Running 1000 simulations.:   0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulation outputs:  tensor([[ 0.6834, -0.2415],\n",
      "        [ 1.3459,  1.5373],\n",
      "        [ 2.1092,  1.9180],\n",
      "        ...,\n",
      "        [ 0.8845,  0.4036],\n",
      "        [ 1.9111,  1.2526],\n",
      "        [ 0.8320,  2.3755]])\n"
     ]
    }
   ],
   "source": [
    "new_theta, new_x = simulate_for_sbi(simulator, restricted_prior, 1000)\n",
    "print(\"Simulation outputs: \", new_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use **all** simulations and run `SNPE` as always:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:root:Found 523 NaN simulations and 0 Inf simulations. They will be excluded from training.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Neural network successfully converged after 118 epochs."
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "13bb699686444c7ca16df8e19a38925c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Drawing 10000 posterior samples:   0%|          | 0/10000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAJgCAYAAABr3Xx6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoFElEQVR4nO3df4zlZ10v8M9nd/bHtNt2ui39wZR2+c0tCOU3uBfFXkmAReBGboIheDFwE40EScQ4YkRDrjq5N5FERYkKVCIqEdGLKYhwFy9aRFsqUCiBtLCUTgtluzu73d3Zn/PcP+bMOp3nO+3Z2TlzZuZ5vZINcz79nvM8Z86w897n+/k+3yylBABAyzYNewIAAMMmEAEAzROIAIDmCUQAQPMEIgCgeQIRANC8kWU8x3X60J4c9gRWyss2/Td/hzVq04UXVrXcunXZrzf7pGuq2ub9h6vamcsv7vs1z1xQ/1o+flk9x4u+ur+v1ytT36tqne/5MTv7er34wYF6jJMnq9rs0aOdT+/6DMrJU3XtVP2aXc9dapzFPj37V4/6d5gVIgCgeQIRANA8gQgAaN5yeogAYE3bfOmlVa2r16VLjtU9P119QJuO1a9XLtje13EREQefVc9x5Hjd4jZybLZjnG1VLY+dqGpnnvWkqrblnrr/aLbj9TbtP1TVTh88WM/v8dfV87u37gtaSm7d0let336h3LK8vjArRABA8wQiAKB5AhEA0DyBCABonqZqANa1riba2SN1A+6mnWNVrVx9eVU7cdloVdv+jXqDw+NPvaqqnR7dXI97unsv0K4G6q2HTle1Y1fWDcYjM3UT9PEnXFSPfbIe48TO8aq27UDdkF2urjdrzI5a6Wga3/SUx1e1iIg4Vb+/rs0eu4xcOlbVZg9OV7XlbrZphQgAaJ5ABAA0zykzWGT35N6Ymp6JiIjxsdG4ZeLGIc8IgEETiGCRqemZ2De5JyIidk3cPOTZALAaBCKaZSWIVnQ1HXfdTXyt6XfH4a5djbvuoB4dta7dmLd17Pjc1UA9/cR6fps6eoZ33NdRjIjZkfoG7Mcvq9/LoSfW3S1Hr6obv7ccqRuox+6uP+cj4/W8tz5Uj/GDZ1xQ1S6+p34vo4eOVbUlHXqoKp3+T9dWtZGv31PVuhqouz7n5TZVC0Q0y0oQAPM0VQMAzROIAIDmCUQAQPP0EAFscOu1gfp85t353K1j9XEXbK9qJx5XH9e1A/WO+85UteM76+NOX9C99tC1g/UPnlU//8Kp+riu5u3DT6xrM1fWO1p3OTFWN2mP7p+tal2N4HnseFXr+r5GdO8M3tVAnRfWDd1dZo/WO5KfObi8nxsrRABA8wQiAKB5AhEA0DyBCABonqZqOEd2uIaV19UE3W+jdddO1V1mD0xXtc1b6udu/9b+qnby2p1V7VhHw3JXo/SRq+pG6aWO7VJedaCqTd9Zz+f0WN3kveNb9a/5o+N1s/TJS+pxz2yv10x23Fsfd+SG8aq29aGOncIjYuRQvQt4l65dqXPHhVVt04V1rXOX8j4IRHCO7HANsPE4ZQYANE8gAgCa55QZG5ZeHwD6JRCxYS3V6zMflMY7dmYF1o6uBuquJtqu3Yo3X3ppn2PUDbjl8rrDeHakPqFyenu9a3NXo/T26bqJOSLiB8+payNH69ecOVE3l29/6qGqduTBenfnI0+pv4fbL6kbm0/fvaOqba43oI6jV9Xfh62Hu99fl03763nH2MV17dhMVTrz/Qf6Hmc5BCKaszAoAUCEQMQGs/g02bzxsdGzq0RWhgBYTCBiQ1lq9Uf/EACPxFVm8AjmV5Z2T+4d9lQAGCArRPAI5leWbMAIsLEJRKxbLquH5ev3thirpd/5dF1R1uXMwYNVresKtU1dVzh95/6qtPmCa6vayPH6th/Hx+qrxE7tqGsREVuue6iqXbNzuqqdOF3/qr7igvq5tx+6rqq99XmfrWp//PXdVa3rqrXDV9X9lpd9rv6cRmbqW4Zsuf9wVYuIKEePVbXZI/Vnej63clkugYh1yy00AFgpAhEbwnyvz2pfQbb46jWrVADrk0DEhjCsILJwXKtUAOuXq8wAgOZZIQJo0Pk0o55vc+tqNHR3jnGyvk1Hl9knXVPVTuzcVh/X8Rt0pON2F6fqu2LMzefOi6ra955ZNygfO1qPve3q01XtZdffWdU+f/AJVe0xFx+pavfdcWVVG32wXjM5fnlVih339b+20tVAnVvr5vRNO8fq5x6Y7nuc5bBCBCvEnkUA65cVIuhDP7f+sGcRwPolEEE8euBx9RjAxiYQQQg8AK0TiAA4J+fbAL3SuxB3PbfvuUzXOypvPlU3LG+94IqqdvLiutn5wPX1GGcuqhulIyJG793cxwwjnn3dd6vaoZP1SvZ/3Xl7Vfutu19Z1S7eVnd+X/y871S1b3xhV1UbfaDedfvwtXWU2P69+nsTEbH5mqur2ulv12NnRwN8V/P1Sjbja6oGAJonEAEAzXPKjDVvo9zEdaO8D4CNSCBizdsoN3HdKO8DYCMSiABY0iB2lV7pBuquZtvO527taL5+zM56Lh3P3XJ/3Xw9clm9bfPOO+um44eu626ePv70mar24+Pfrmqf+ebTqtrs8fo1f+t03UD9zid+oqr9/r316vSJM3Uc2PX8e6vaPV+od/G+7M66CX3T/kNVLSKiXFJv293vZ9rvTuPLJRDBAC0+TQbA2iQQwQpbvMnj/GkyANYugQhWmGZpgPXHZfcAQPOsEAGwpPPZVXqpHaRXcnfhiIjZo0er2uZLL60P7Lf5+li9k/Opa+sG6pM76jWF09s7Xq/uOY6IiJHv1Af/w/GnV7U3v+Cfq9oHPvvSqvba8S9XtXfe+dqq9iPjd1e1p43eX9X+17+9vKpte8qRqnZf1I3S187UO3tHRBy/rP6ZuGj/WOexi5UD030dt9ydywUi1pV+7joPAOdKIGJd2Wj9OTZrBFgbBCIYIps1AqwNmqoBgOZZIQLgnPTbFL3Ucf02vXY9v6vW1UA9e6Sj0fqaq6vamW9/t6/jTl1cz/ni79TN18eu3FbVjo7Xu1dHROz4Tl07OF7vk/3RP6lPpZdn1t+Hbxy7sqo9/6p7qtoFm+rn/vuRa6va//7hv6pqv/Tpn6pqmzs24r7vh7t7PB/3mYfq4iUXVaUy9b26tsLN+ItZIQIAmicQAQDNE4gAgOYJRABA8zRVA6xT/e4YvdbG6HdX600XXtjX6505eLCv1yvTh+sxdo7Vx22pfzVue7BuoP7ei+tm4K5dqbceqmsR3c3W5XRdO3LtbFV7+pOmqtqzd9QN1Ndvr4/78P4XV7UHT9Tf67/Z/5x6fpvrpu/tB+o5H73uTFWLiJi5sm623rG//lzyinpn8HigLpWTp+raMn8+BSJYIxbvwm2TRoDVIxDBGrEwANmkEWB16SECAJonEAEAzXPKDGCdGvTOvec7xlJN0bNH612kc+uWeuyOhtl+x8mtHbthd4yRW+panKo7ozcdOlbVHvv3de3BF11Rv96RuhQRceH36gblw8f728X7zn2PrWrjF0xXtT+66z9XtSdc+mBV+6kr/62qfedk3dj8Z6/+x6r2lP/336taOdrxfY2IY1fU9ZGZepxtt3+rqm3qaLQux2bqWscu5f2wQgQANE8gAgCa55QZDMH8JfbjY903QHQJPsDqEohgCB4t4LgEH2B1CUQADF1XE3Tu6GjK7rPRurOB+oJ6RbZr9+p4bN0YffqS7X0NO3K8bpQeOVbvNB0RcXzn5qq2qWOD565dn3d8ZVtV+7/7n1XVLr9+f1X74jd2VbUbLrm3ql2wqW6o/5/7n1bVzpzpv/vm+OX1rtYjx+vPavaFT6pq2//pzqrW9TPSbzP+YnqIAIDmCUQAQPMEIgCgeQIRANA8TdWsWbsn98bU9MySl6YDtdxSNyevxo7WXeN27Ui9lNmO3YW75t01zqadY/Vzu3YrvuSi+vUuvKCqnepooB45dLw+7vL6uWO33l/P5YKlGrIvqSond9S/lrdM12sXDz2l3k370sceql/vdN24veXCuun4FRd9papNfOsnq9rTx+r39+InfLuqff7fn1rVIiJm6+nE0avq97fzc9+tamXXNfXrfbMee7kEItasqemZ2De5Z9jTAKABTpkBAM0TiACA5glEAEDz9BABbCCr0UA9iHG7nr/pwv52Ie5qoO7awThP1Y3IXQ3PW+6pd3eevbxugN56z4Gqdubyi6vaicu6m6oPX1v/Cj69oz6uqxF59zO/WdW++pfX12M8sd4l+/U/8vmq9j/ueGNV27Gtv8/0zv1XVrWyrWPL7YiIqN/Mtul6d+/YUn9v8ljd2J4dO5LbqRoAYJkEIgCgeQIRANA8gQjWuPGx0dg1cXPsntw77KkAbFiaqmGNu2XixoiY27l718TNETEXkubrsN507Tbdpd/m2K4dsTd3NFV3NVB3NerOHpyuB+mode6cvKVeZxj9Vt18PWdnVTl5cd10fPKSrGpf+7O6gfrEjz1U1S4cqZubj83W3/9TZzo6tzvsP1F3fXc1X89+pfszPtPRXz52d/0ZdO7u3dEUv+nSsfq5p5bXVC0QwTqxMADNByMAVoZTZgBA8wQiAKB5AhEA0DyBCABonqZqANaNfm/xMXtguqplxy0+znRdoXbppfVzL7ygnst9D1S1kbiiqp26ur6dR0TEBV+9ry4+47FVaeuh+vYUXVdrbbr9oqp29Kr61h37r6yvFDt6X/3c//Ki+vYg/+fW51S1bd+vo8SJ6+srwiIirvmH+oq52ZGOtZmO72103aajY4yuz74fVogAgOZZIYJ1aH6zxvmv7UkEcH4EIliH7EkEsLKcMgMAmmeFCIBV1dUY3aXrFh+dta11reu2H13HZcdxZw4erI/raMje1HF7kK5bgWy9+0hVi4g4c3V9647R7x6uajtHLqlqB66vf33vmKpbjLceqtc9bpuub/vRdaONT/7986vapT9U34bk9Dcvq4/7Rvd6y6aTdZP31gdnqlqO1Y3oZbr+3nQ1WmdHrR9WiACA5glEAEDzBCIAoHkCEQDQPE3VAAxEVwN0RHdTddex/R7X1QTd5czB/pq5u3a+7mrSnu1qtO56wcfWu1dHRGw6dKyv+Wx9qB57/LN1I/Lslnr0B59R77B9yd31GNsPnOlrLsfurhvBtx2vG6VHjnftIR0x+v2OBur791e109/v2Km6Q3eTvaZqAIBlEYgAgOYJRABA8wQiAKB5mqoBOG/9NkWfr/PZ5brzuI4G3K4G6i5dO1WXk/X8Zr/2jc7nb76yu9l6sdNPqBuZNx87XdVGDtW7ZF/xb/V7OX5VPe/t36sbxA8/pd4t+rI76l23O5vDf1DvaB0RkRfWTd5d3+/Nl15a1bqa2Ls+v9mj9XH9sEIEADTPChFryu7JvTE1PXdZ5vjY6JBnA0ArBCLWlKnpmdg3uWfY0wCgMU6ZAQDNs0IEwHk73wbqlW7A7vf1+j2ue/fqjp20OxqtN3fUIiKio5k4x+pG5tF//049dkeDcXSMk5fXzcldr5db6ubksUP1GOVYvdN01/so111dHxcR2fGa0dUYfWC68/nVcR0N1P021C9mhQgAaJ5ABAA0TyACAJonEAEAzdNUDevc+Nho7Jq4+ezXt0zcOOQZwdrU1Wzb707H/Tbvlo5m4K4xljJ77/19Pb9rd+euscv3H6hfr8+m466duLt07s79pTu7j+1zR/Pz2fl8uQ36AhGscwsD0HwwAuDcOGUGADTPChFD53YdAAybQMTQuV0HAMMmEAHQhK5m29XYIXuppup+G7X7Pa6z+brP99fduN3x/epq5u6zKfpcj10uO1UDACyTQAQANE8gAgCaJxABAM3TVA0b0OKtDOxeDcPT1RR9vs6nQbzf5w5it+jzec1+Lff1rBDBBrLwNh77JvfEvsk9Z4MRAEuzQgQbiJUggOWxQgQANE8gAgCa55QZACzDSjcDn+9rrnTD8lp7f4NmhQgAaJ5ABAA0TyACAJonEAEAzdNUDQAbwFpuWF4PBCKGYvGtJQBgmAQihmJqeib2Te4Z9jQAICIEItjwFt7fzI1eAboJRLDBLQxA88EIgIdzlRkA0DyBCABonkAEADRPIAIAmicQAQDNE4gAgOYJRABA8wQiAKB5AhEA0DyBCABonkAEADRPIAIAmicQAQDNE4gAgOYJRABA8wQiAKB5AhEA0DyBCABonkAEADRPIAIAmicQAQDNE4gAgOYJRABA80aGPQFg9YyPjcauiZvPfn3LxI1DnhHA2iAQQUMWBqD5YASAU2YAAFaIWF27J/fG1PRMjI+NDnsqAHCWQMSqmpqeiX2Te4Y9DQB4GKfMAIDmCUQAQPMEIgCgeQIRANA8gQgAaJ5ABAA0z2X30Ci38QD4DwIRNMptPAD+g1NmAEDzBCIAoHkCEQDQPIEIAGieQAQANE8gAgCaJxABAM0TiACA5tmYkYHbPbk3pqZnImJuR2QAWGsEIgZuanom9k3uGfY0AGBJTpkBAM0TiACA5glEAEDzBCIAoHkCERDjY6Oxa+Lm2D25d9hTARgKV5kBccvEjRERsWvi5iHPBGA4BCLgrPmVovmv54MSwEYnEAFnLQxAVouAlughAgCaJxABAM0TiACA5glEAEDzBCIAoHkCEQDQPIEIAGiefYgYiN2Te2NqeiYi5jb4A4C1TCBiIKamZ2Lf5J5hTwMA+iIQsWKsCgGwXglErBirQgCsV5qqAYDmCURAp/GxuTvfu8kr0AKnzIBOC+98D7DRWSECAJonEAEAzROIAIDmCUQAQPMEIgCgeQIRANA8gQgAaJ5ABAA0TyACAJonEAEAzROIAIDmCUQAQPMEIgCgeQIRANA8gQgAaJ5ABAA0TyACAJonEAEAzRsZ9gRY/3ZP7o2p6ZkYHxsd9lQAYFkEIs7b1PRM7JvcM+xpAMCyOWUGADRPIAIAmicQAQDNE4gAgOYJRABA8wQiAKB5Lrunb/P7DUVEjI+Nxi0TNw55RgCwMgQi+rZwv6FdEzcPeTYAsHKcMgMAmmeFiGUZHxs9u0rklh0ArHcCEcuifwiAjUQg4hEtbqQGgI1IIOIRuXErAC3QVA0ANE8gAgCa55QZFX1DALRGIKKibwiA1jhlBgA0zwoRZ82fKnOaDIDWZCll2HMAABgqp8wAgOYJRABA8wQiAKB5AhEA0LxzvsosM78aEccHMJd+XR4R+xsdv+X3PuzxW37vERHbSynPGOL4AAO1nMvuj5dSnrfiM+lTZt7W6vgtv/dhj9/ye58ff1hjA6wGp8wAgOYJRABA85YTiP5oxWdh/PUwduvjt/ze18L4AANlp2oAoHlOmQEAzROIAIDmnVMgysw3ZOZXMvOOzPx8Zj5rUBPLzKdl5r9k5onMfMcjHHdTZn47M7/U+3PDMOezQmNlZv5uZt7V+34/Z4nj/jEzv7HgvV8xwDm9vDfWXZk5MaxxMvNNmfmDBe/5LYOaS2+8D2TmA739twbq0cbKzJdm5qEF7/1dA57P4zLzs5l5Z2Z+LTN/YZDjAQzTue5D9O2I+NFSysHMfEXMNVq+cOWnFRERByLibRHx2j6O/aVSykcHNI955zKf8/WKiHhy788LI+IPY+nv8xtKKQPdIyYzN0fEeyPiZRFxb0TcmpkfL6XcOaRxPlJKeetKjv0IboqI34+ID62Rsf6plPKqVZhLRMTpiPjFUsrtmXlRRHwxMz+90p87wFpwTitEpZTPl1IO9h5+ISKuWfkpnR3rgVLKrRFxalBjnItVns9rIuJDZc4XImIsM69ehXGX8oKIuKuU8q1SysmI+MveHNfrOH0rpXwu5sLwhhqrH6WU+0spt/e+figivh4R48OdFcBgnE8P0Zsj4pMrNZHz9Ju9U0vvycxtw57MChiPiO8ueHxvLP2L6IO90ye/lpm5BuazGuP8ZO/z/mhmPm4A81jLXpyZX87MT2bm01dr0MzcFRHPjoh/Xa0xAVbTsgJRZv5YzAWiX17Z6SzLr0TE0yLi+RGxM9bGnFbLG0opPxQRL+n9eeOQ57Ma/i4idpVSnhkRn46IPx3yfFbT7RFxXSnlWRHxexHxt6sxaGbuiIi/joi3l1IOr8aYAKvtUQNRZv78gibOx2bmMyPiTyLiNaWUB1dyMovH6uc5vWX9Uko5EREfjLnTLkObz0qMFRH3R8TClY9rImJq8XNKKVO9/30oIv48VvC9LzLVz3xWY5xSyoO9zzpi7ufwuQOYx5pUSjlcSjnS+/oTEbElMy8f5JiZuSXmwtCHSykfG+RYAMP0qIGolPLeUsoNpZQbYq4J+2MR8cZSyjdXejILxyql3NfPc+Z7a3qni14bESt2NdBy5rMSY8Xcv/x/une12Ysi4lAp5f6Fx2fmyPwvw94vrVfFCr73RW6NiCdn5uMzc2tEvD4iPj6McRb1Ur065vpampCZV82fFs3MF8Tc/39X9B8li8bLiHh/RHy9lPI7gxoHYC0416vM3hURl0XEH/T+Xj49qDtwZ+ZVEXFbRFwcEbOZ+faIuL6UcjgzPxERb+mFlA9n5mMiIiPiSxHxs6s9nwEM94mIeGVE3BURxyLiZxbM40u90LQtIj7VC0ObI+IzEfHHA5hLlFJOZ+ZbI+JTvbE+UEr52mqNk5nvjojbSikfj4i3ZearY+4KqAMR8aaVnsdCmfkXEfHSiLg8M++NiF8vpbx/tcaKiC0REaWU90XE6yLi5zLzdETMRMTry2C3mt8dc6dh7+itXEZEvLO3OgWwobh1BwDQPDtVAwDNE4gAgOYJRABA8wQiAKB5AhEA0DyBiMjM38jMd/S+fndm/vh5vNaq3R0eAFaKQMTDlFLeVUr5zHm8xE0R8fIVmg4ArAqBqFGZ+auZ+c3M/OeIeOqC+k2Z+bre1/sy87d7txS5LTOfk5mfysy7M7NzA8y1dsd2AOjHue5UzQaQmc+Nudti3BBzPwO3R8QXlzj8nlLKDZn5nphb/dkdEdtj7jYh7xv4ZAFgFQhEbXpJRPxNKeVYRERmPtJ9yeb/2x0RsaN3I9mHMvNEZo6VUqYHO1UAGDynzHg083eWn13w9fxjgRqADUEgatPnIuK1mTmamRdFxE8Me0IAMEwCUYNKKbdHxEci4ssR8cmIuHWlXrt3x/Z/iYinZua9mfnmlXptABgUd7sHAJpnhQgAaJ5ABAA0TyACAJonEAEAzROIAIDmCUQAQPMEIgCgeQIRANA8gQgAaJ5ABAA0TyACAJonEAEAzROIAIDmCUQAQPMEIgCgeSPDngDAKivDngCwqrKfg6wQAQDNE4gAgOYJRABA8wQiAKB5AhEA0DyBCABonkAEADRPIAIAmicQAQDNE4gAgOYJRABA89zLDAAatntyb0xNz5x9PD42GrdM3DjEGQ2HQAQADZuanol9k3vOPt41cfMQZzM8AhEAbGCLV4AWGx8bXcXZrF0CEQBsYItXgB7N+Njow1aJWjmFJhABAGctDj+tnEITiABgA+lqkubRCUQAsIGc6yky5tiHCABonkAEADRPIAIAmicQAQDNE4gAgOYJRABA8wQiAKB5AhEA0DyBCABonkAEADRPIAIAmicQAQDNE4gAgOYJRABA80aGPQEAYPl2T+6NqemZs4/Hx0aHOJv1SyACgHVsanom9k3uGdjrj4+Nxq6Jmx/2+JaJGwc23rAIRADAkhaHn4XhaCPRQwQANE8gAgCaJxABAM0TiACA5glEAEDzBCIAoHkCEQDQPIEIAGieQAQANE8gAgCaJxABAM0TiACA5rm5KwCsYbsn98bU9MzZxxv1bvPDJhABwBo2NT0T+yb3nH28Ue82P2xOmQEAzROIAIDmCUQAQPP0EAHAOjI+NvqwPqLxsdEhzmbjEIgAYB0Z9hVmXYFs2HNaCQIRANC3xeFno1z1pocIAGieQAQANE8gAgCaJxABAM0TiACA5glEAEDzBCIAoHkCEQDQPIEIAGieQAQANE8gAgCaJxABAM0TiACA5glEAEDzBCIAoHkCEQDQPIEIAGieQAQANG9k2BMAAP7D7sm9MTU9c/bx+NjoEGfTDoEIANaQqemZ2De5Z9jTaI5TZgBA8wQiAKB5AhEA0DyBCABonkAEADRPIAIAmicQAQDNsw8RAAyRjRjXBoEIAIbIRoxrg1NmAEDzBCIAoHkCEQDQPIEIAGieQAQANE8gAgCaJxABAM0TiACA5glEAEDz7FQNACzb+Nho7Jq4+WGPb5m4cYgzWh6BCABYtsXhZ2E4Wk+cMgMAmicQAQDNE4gAgOYJRABA8wQiAKB5AhEA0DyBCABonkAEADRPIAIAmicQAQDNE4gAgOYJRABA8wQiAKB5AhEA0DyBCABonkAEADRPIAIAmicQAQDNE4gAgOYJRABA8wQiAKB5AhEA0DyBCABonkAEADRPIAIAmicQAQDNE4gAgOYJRABA8wQiAKB5AhEA0DyBCABonkAEADRvZNgTAICW7J7cG1PTM2cfj4+NDnE2zBOIAGAVTU3PxL7JPcOeBos4ZQYANE8gAgCa55QZALBixsdGY9fEzQ97fMvEjUOcUX8EIgAYoNaaqBeHn4XhaC0TiABggDRRrw96iACA5glEAEDzBCIAoHkCEQDQPIEIAGieQAQANE8gAgCaJxABAM0TiACA5glEAEDzBCIAoHkCEQDQPIEIAGieu90DwAraPbk3pqZnzj4eHxsd4myGb3xsNHZN3Pywx7dM3DjEGXUTiABgBU1Nz8S+yT3DnsaasTj8LAxHa4lTZgBA8wQiAKB5AhEA0DyBCABonkAEADRPIAIAmicQAQDNE4gAgObZmBEAzoOdqTcGgQgAzoOdqTcGp8wAgOYJRABA8wQiAKB5AhEA0DyBCABonkAEADRPIAIAmicQAQDNE4gAgOYJRABA8wQiAKB5AhEA0DyBCABonkAEADRPIAIAmicQAQDNGxn2BABgPdk9uTempmfOPh4fGx3ibFgpAhEAnIOp6ZnYN7ln2NNghTllBgA0TyACAJonEAEAzdNDBACsmvGx0dg1cfPDHt8yceMQZzRHIAIAVs3i8LMwHA2TU2YAQPMEIgCgeQIRANA8gQgAaJ5ABAA0TyACAJonEAEAzROIAIDmCUQAQPMEIgCgeQIRANA8gQgAaJ5ABAA0TyACAJonEAEAzROIAIDmCUQAQPNGhj0BAFjLdk/ujanpmbOPx8dGhzgbBkUgAoBHMDU9E/sm9wx7GgyYU2YAQPOsEAEAQzM+Nhq7Jm5+2ONbJm5c9XkIRADA0CwOPwvD0WpyygwAaJ5ABAA0TyACAJonEAEAzROIAIDmCUQAQPMEIgCgefYhAoAF3LusTQIRACzg3mVtcsoMAGieQAQANE8gAgCaJxABAM0TiACA5rnKDICmucyeCIEIgMa5zJ4Ip8wAAAQiAACBCABonkAEADRPUzUAsGaMj43GrombH/b4lokbBz6uQARAU1xmv7YtDj8Lw9EgCUQANMVl9nTRQwQANE8gAgCa55QZALBmrVaTtUAEAKxZq9Vk7ZQZANA8gQgAaJ5ABAA0TyACAJonEAEAzROIAIDmCUQAQPPsQwTAhuZmrhvLoDZqFIgA2NDczHVjGdRGjU6ZAQDNE4gAgOY5ZQbAhqJniOUQiADYUPQMsRwCEQDrmhUhVoJABMC6ZkWobYsvw1+s358NgQgAWLdWYg+iCIEIgHXGKTIGQSACYKi6As4j/avfKTIGQSACYKgWB5zdk3sfsSfEihCDIBABMFDnugK0Uj0hcC4EIgAG6tFWgKz4sBYIRADr3OIVmEezeIXmXJ9/rhYHHitArEVZShn2HAAAhsrNXQGA5glEAEDzBCIAoHkCEQDQPFeZAU3JzK9GxPEhDX95ROwf0titj9/yex/2+MN+79tLKc94tIMEIqA1x0spzxvGwJl527DGbn38lt/7sMdfC++9n+OcMgMAmicQAQDNE4iA1vxRo2O3Pn7L733Y46+L926nagCgeVaIAIDmCURAEzLzDZn5lcy8IzM/n5nPGuBYT8vMf8nME5n5jkc47qbM/HZmfqn354ZhzmeFxsrM/N3MvKv3/X7OEsf9Y2Z+Y8F7v2KAc3p5b6y7MnNiUOP0M1Zmvikzf7Dgfb9lwPP5QGY+0NtuYqAebazMfGlmHlrw3t814Pk8LjM/m5l3ZubXMvMXHul4l90Drfh2RPxoKeVgZr4i5voKXjigsQ5ExNsi4rV9HPtLpZSPDmge885lPufrFRHx5N6fF0bEH8bS3+c3lFL6uiR6uTJzc0S8NyJeFhH3RsStmfnxUsqdQxzrI6WUt670+Eu4KSJ+PyI+tEbG+qdSyqtWYS4REacj4hdLKbdn5kUR8cXM/PRSn70VIqAJpZTPl1IO9h5+ISKuGeBYD5RSbo2IU4Ma41ys8nxeExEfKnO+EBFjmXn1Koy7lBdExF2llG+VUk5GxF/25rjex+pLKeVzMReIN9RY/Sil3F9Kub339UMR8fWIGF/qeIEIaNGbI+KTw55Ez2/2Ti29JzO3DXsyK2A8Ir674PG9sfQvoQ/2Tp38WmbmGpjPao31k73P/KOZ+bgBzWWtenFmfjkzP5mZT1+tQTNzV0Q8OyL+daljBCKgKZn5YzEXiH552HOJiF+JiKdFxPMjYmesjTmtljeUUn4oIl7S+/PGIc9ntfxdROwqpTwzIj4dEX865Pmsptsj4rpSyrMi4vci4m9XY9DM3BERfx0Rby+lHF7qOIEI2LAy8+cXNHA+NjOfGRF/EhGvKaU8OMix+nlOb0m/lFJORMQHY+6Uy9DmsxJjRcT9EbFw1eOaiJha/JxSylTvfx+KiD+PFXzvi0z1M5/VGquU8mDv846Y+1l87oDmsuaUUg6XUo70vv5ERGzJzMsHOWZmbom5MPThUsrHHulYgQjYsEop7y2l3FBKuSHmLiL5WES8sZTyzUGOVUq5r5/nzPfW9E4XvTYiVuxKoOXMZyXGirl/9f9072qzF0XEoVLK/QuPz8yR+V+EvV9Yr4oVfO+L3BoRT87Mx2fm1oh4fUR8fFhjLeqnenXM9bU0ITOvmj81mpkviLkMsqL/MFk0XkbE+yPi66WU33m0411lBrTiXRFxWUT8Qe/v5NODuuFkZl4VEbdFxMURMZuZb4+I60sphzPzExHxll5I+XBmPiYiMiK+FBE/u9rzGcBwn4iIV0bEXRFxLCJ+ZsE8vtQLTdsi4lO9MLQ5Ij4TEX88gLlEKeV0Zr41Ij7VG+sDpZSvreZYmfnuiLitlPLxiHhbZr465q6AOhARbxrEXOZl5l9ExEsj4vLMvDcifr2U8v7VGisitkRElFLeFxGvi4ify8zTETETEa8vg90denfMnYq9o7d6GRHxzt7qVD1/O1UDAK1zygwAaJ5ABAA0TyACAJonEAEAzROIAIDmCUQArFuZ+RuZ+Y7e1+/OzB9f5uuc053R2XjsQwTAhlBKedd5PP2c7ozOxmOFCIB1JTN/NTO/mZn/HBFPXVC/KTNf1/t6X2b+du+WIrdl5nMy81OZeXdmVhtgnuud0dl4rBABsG5k5nNj7pYYN8Tc77DbI+KLSxx+Tynlhsx8T0TcFHM7F2+PuduEvO8RxtgVj3JndDYegQiA9eQlEfE3pZRjERGZ+Uj3JZv/b3dExI7eys9DmXkiM8dKKdOLn9DvndHZeJwyA2Cjmr+r/OyCr+cfVwsC53JndDYegQiA9eRzEfHazBztNT//xEq86LneGZ2NRyACYN3oNT5/JCK+HBGfjIhbV+il5++MfmOvEftLmfnKFXpt1gF3uwcAmmeFCABonkAEADRPIAIAmicQAQDNE4gAgOYJRABA8wQiAKB5AhEA0Lz/DzPKZgMwsMoIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "restriction_estimator.append_simulations(\n",
    "    new_theta, new_x\n",
    ")  # Gather the new simulations in the `restriction_estimator`.\n",
    "(\n",
    "    all_theta,\n",
    "    all_x,\n",
    "    _,\n",
    ") = restriction_estimator.get_simulations()  # Get all simulations run so far.\n",
    "\n",
    "inference = SNPE(prior=prior)\n",
    "density_estimator = inference.append_simulations(all_theta, all_x).train()\n",
    "posterior = inference.build_posterior()\n",
    "\n",
    "posterior_samples = posterior.sample((10_000,), x=torch.ones(2))\n",
    "_ = pairplot(posterior_samples, limits=[[-2, 2], [-2, 2]], fig_size=(3, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further options for tuning the algorithm\n",
    "\n",
    "- the whole procedure can be repeated many times (see the loop shown in \"Main syntax\" in this tutorial)  \n",
    "- the classifier is trained to be relatively conservative, i.e. it will try to be very sure that a specific parameter set can indeed not produce `valid` simulation outputs. If you are ok with the restricted prior potentially ignoring a small fraction of parameter sets that might have produced `valid` data, you can use `restriction_estimator.restrict_prior(allowed_false_negatives=...)`. The argument `allowed_false_negatives` sets the fraction of potentially ignored parameter sets. A higher value will lead to more `valid` simulations.  \n",
    "- By default, the algorithm considers simulations that have at least one `NaN` of `inf` as `invalid`. You can specify custom criterions with `RestrictionEstimator(decision_criterion=...)`"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
